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We show that discrete quasiprobability distributions defined via the discrete Heisenberg-Weyl 
group can be obtained as discretizations of the continuous SU{N) quasiprobability distributions. 
This is done by identifying the phase-point operators with the continuous quantisation kernels 
evaluated at special points of the phase space. As an application we discuss the positive-P function 
and show that its discretization can be used to treat the problem of diverging trajectories. We study 
the dissipative long-range transverse-field Ising chain and show that the long-time dynamics of local 
observables is well described by a semiclassical approximation of the interactions. 
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I. INTRODUCTION 

Phase-space representations of quantum mechanics, 
usually called quasiprobability distributions, provide a 
natural language for the quantum-classical correspon¬ 
dence of non-relativistic quantum mechanics. They were 
first developed in the context of harmonic oscillators [1] 
and extensively applied in quantum optics 013 ]. It has 
been realised that the quantum-classical correspondence 
is intimately related to symmetry properties of the un¬ 
derlying classical phase-space. This led to an axiomatic 
approach to the phase-space formulation of quantum me¬ 
chanics known as Stratonovich-Moyal-Weyl correspon¬ 
dence 00. In order to apply the ideas of Stratonovich 
and Berezin to finite dimensional Hilbert spaces (e.g. 
spin systems) two main research directions are followed. 
First is to construct a continuous phase space with a 
natural continuous symmetry of the system (e.g. SU(2) 
for spin systems) jbUTQ]. In many cases this permits to 
rewrite the evolution equations in terms of stochastic dif¬ 
ferential equations mm or even ordinary differential 
equations [HIS] that can be efficiently simulated. The 
second direction is to exploit the discrete nature of the 
system and to define phase-space distributions on the ba¬ 
sis of the discrete Heisenberg-Weyl group [HHii]. This 
formulation enables a simpler representation of quan¬ 
tum states and has been applied in quantum information 
[23l l24] and quantum tomography [25l |26] . 

We show that the two approaches to formulate phase- 
space quasiprobability distributions on finite dimensional 
Hilbert spaces are equivalent, in the sense that discrete 
distributions can be obtained by evaluating the continu¬ 
ous ones at special points of the phase space. This iden¬ 
tification provides a formal justification of Monte-Carlo 
methods on finite dimensional phase spaces proposed in 
[27l l28] , their systematic expansion beyond the semiclas¬ 
sical (truncated Wigner) approximation [T2j |29] and ex¬ 
tension to open quantum systems. Hence, the revealed 
relation between continuous and discrete phase spaces 
should be relevant to further development of contempo¬ 
rary phase-space methods for simulation of long-range 


many-body quantum systems. A similar discretisation 
was obtained in case of the SU (2) group [30l [3T] by con¬ 
tracting a continuous SU(2) kernel to a discrete one. Here 
we take the opposite route and extend the discrete kernel 
to the continuous one for any dimension of the Hilbert 
space N. 

In order to prove the main result of the paper we first 
review the basic properties of the discrete and SU{N) 
quasiprobability distributions on A'-dimensional Hilbert 
space. 


A. Discrete quasiprobability distributions 

A one parameter family of discrete Weyl symbols for 
an operator A acting on a A'-dimensional Hilbert space 
is usually defined as 

trW*k (1) 

with the parameter s denoting the ordering of the dis¬ 
tribution; s = 1,0,—1 for normal, symmetric, and anti¬ 
normal ordering, respectively. The phase-point operators 
/3’ /^ = 0,1, 2,... — 1 satisfy the following prop¬ 

erties usmuEi]: 

la) Hermicity; 

2a) Normalization: tr = 1, 

3a) Covariance: A^*!^ where 

denotes a unitary irrep of the discrete Heisenberg- 
Weyl group of order N, 

4a) Traciality: tr A^^^gA^^T^l = N 6 a,a'Si 3 ,/ 3 '- 

The quasiprobability distribution is given by the Weyl 
symbol of the density matrix. For each satisfy¬ 

ing the properties la-4a the Weyl symbol 0 defines a 
one-to-one mapping between operators on a finite dimen¬ 
sional Hilbert space and functions on a discrete phase 
space. The explicit form of the phase-point operators is 
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not relevant for our purpose and will be omitted; we re¬ 
fer the interested reader to references [m [151 [m (foi* a 
particular example see Section Ilk. The kernel can 

be rewritten in the following us^ul form m 

( 2 ) 

The decomposition (§ is valid for any unitary realisation 
of the discrete Heisenberg-Weyl group of order N. 

The constant matrix Aq^q is determined by demanding 
that for 5 = —1 the quasiprobability distribution repre¬ 
sents the discrete Q function. This gives the boundary 
condition 

=Ta,0\O){O\Tl^, (3) 

where |0) denotes a chosen vacuum state. Other Aq^q 
can be determined by using the traciality condition [21] . 
Phase-point operators satisfying the properties la-4a can 
be defined and calculated for any Hilbert space dimen¬ 
sion N. If in addition N is a prime number or a power 
of a prime number the quasiprobability distribution is 
defined on a discrete phase space with a well defined ge¬ 
ometry [laiiHiiia. For such systems (sometimes called 
Galois quantum systems) phase-space methods based on 
discrete symplectic transformations were developed that 
are similar to the ones used for harmonic oscillator sys¬ 
tems m- Phase spaces with non-prime number dimen¬ 
sions N do not have a well defined geometry. As a con¬ 
sequence the standard construction of the mutually un¬ 
biased basis does not work, leaving this fundamental and 
interesting problem unsolved in the Hilbert spaces of non¬ 
prime number dimension. In addition, on even dimen¬ 
sional phase spaces the discrete kernels are not uniquely 
determined OEl]. This, however, does not affect our 
results. 


B. SU{N) quasiprobability distributions 

The SU{N) quasiprobability distributions on a mani¬ 
fold M. are constructed in a similar manner as the dis¬ 
crete quasiprobability distributions through the quanti¬ 
sation kernel The Weyl symbol is given by 

= tvAA^^\n), (4) 

with the parameter s = —1,0,1 denoting the normal, 
symmetric, and anti-normal ordering. The quantisation 
kernel has to satisfy similar properties as in the discrete 
case, namely [4] 

Ib) Hermicity: (A^")(fl))t = 

2b) Normalization: tr A^^^(fl) = I, 

3b) Covariance: = A{g)A^^\Q)A'^{g), where 

A{g) denotes a unitary irrep of a coset element 
geg = SU{N)/U{N-l), 


4b) Traciality: d/i(fl)tr (A*^^^(fI)A*^ = 

where d/i(fl) is the invariant measure, 
denotes a differentiable 5-ordered function 
on the manifold Ai being isomorphic to g. 

For each quantisation kernel satisfying the prop¬ 

erties Ib-4b the Weyl symbol Q represents a one-to- 
one mapping between operators on the Wdimensional 
Hilbert space and smooth functions on the classical man¬ 
ifold M. There are several constructions of the quantisa¬ 
tion kernel for the SU{N) group [6l[9j[T0]. Here we focus 
on the fundamental representation for which the most 
explicit expression for the quantisation kernel was given 
in uni- The explicit form of the SU{N) kernel shall be 
omitted; interested reader is referred to [6l [9] [10] (for a 
particular example see Section [Ilj). As in the discrete case 
all realisations of the SU{N) group admit the following 
decomposition of the kernel [9] 

(5) 

where is a constant diagonal matrix containing the 
essential information about the quasiprobability distribu¬ 
tion. It can be determined by demanding that for s = — 1 
the usual Q function should be recovered, namely 

= A(fi)|0)(0|A(^^)r (6) 

All other are determined by the traciality condition 

mm- 


II. CONSTRUCTION OF THE SU{N) KERNEL 
FROM THE PHASE-POINT OPERATORS 

In the following we shall prove that continuous SU(N) 
quasiprobability distributions for the fundamental repre¬ 
sentation of the SU{N) group can be constructed from 
the discrete phase-point operators. We start by realising 
that the discrete Heisenberg-Weyl group of order N is 
a subgroup of the coset g. Hence, the realisations 
associated to the Wdimensional Hilbert space can be 
seen as fundamental realisations of particular group ele¬ 
ments ^ namely = A(flQ,^^), suggesting that 
the phase-point operators can be regarded as the SU(N) 
kernels evaluated at phase-space points This is im¬ 

mediately clear for the Q function, however, it can be 
shown for any s without knowing the precise form of the 
kernels A^^^^ and A^^^(fl). 

The discrete phase-point operators can be extended to 
the SU{N) group by defining 

A«yn) = A(lI)Ai>t(^). (7) 

The extended Weyl symbol for an operator A is then 
defined as 

w[^\a,l3,n)=tT[AAl:^^{n)). ( 8 ) 
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In the following we shall show that for any fixed Q the ma¬ 
trix (a, /3, Q) represents a discrete Weyl symbol and 

that for any fixed pair a, [3 the function (a, rep¬ 
resents a continuous Weyl symbol. The properties la, 2a, 
and 4a of the phase-point operators are trivially satisfied. 
The covariance property 3a follows by using a rotated re¬ 
alisation of the discrete Heisenberg-Weyl group, namely 

Further, the properties Ib- 
3b of the continuous quantisation kernel are evident. Fi¬ 
nally, the continuous traciality property 4b follows from 
the invariance of the Haar measure and from the fact that 
Aa,/3(^) are also discrete phase-point operators for any 
but fixed ft e Q. By using these two properties any op¬ 
erator A acting on the Hilbert space can be decomposed 
as 

A = (9) 

a, 13 

= 1 /■ tr 

Jm ^ 

Jm 


where ft^ denotes the invariant volume of the phase 
space M. From above equality the continuous tra¬ 
ciality condition follows 

( 10 ) 

Hence, the kernel AcK^/3(f2) can be used to generate a 
discrete or a continuous phase space. 

As an example we consider a spin 1/2 case, where the 
phase space is a sphere, which can be useful to gain some 
intuition on the mapping between the discrete to the con¬ 
tinuous phase space. We start by choosing the following 
phase-point operators 


+ 0 \ 

An , — 1 

V 0 pi - ^/3) ) 

5 ^0,1 — 1 


A-1,0 — 


Ao,0 = 


1(3-V3) (-1; V3^ 

H3 + \/3) 




1(3-V3) ^ 3 . 

i(3 + V3) 


( 11 ) 


Ai.i = 


1 ( 3 - 73 ) - 7 = 1 , 




(_l)2/3 /2^-i^ 1(3 + 73) 


Using the SU{2) kernel for the Wigner function 

/ v^(i-hP) , 1 V3z 

2{\z\^ + l) ^ 2 |d" + l 

Vsz v^(hl''-i) , 1 

V hl"+i 2(|d"+i) 2 

one can easily check that the phase-point operators 
© correspond to the continuous kernel © evaluated 
at special points of the phase space, namely Aij = 
A(°3zi,j) with 2:0,0 = 0, 20,1 = \/2expiip, 21,0 = 
(_l)2/372expi+ and 2:1,1 = —(—l)^A72expi+. Apply- 
ing the stereographic projection 2) = tan(^/2) exp(i(/)) we 
find that the points corresponding to the phase-point op¬ 
erators form vertices of a regular tetrahedron embedded 
in the continuous SU{2) phase space (see Figure [^. In 
Appendixwe use a slightly more general discrete phase 
space points, namely we also allow the rotation around y 
axis of the phase space. 

We have established a connection between discrete 
quasiprobability distributions and SU{N) quasiproba¬ 
bility distributions on finite-dimensional Hilbert spaces. 





Figure 1 : Points in the continuous SU{ 2 ) phase space cor¬ 
responding to a particular choice of phase-point operators 
©■ The red, green, blue and black points correspond to 
2^0,0, 2^0,1, 2^1,0 and zi,i, respectively. Changing (p rotates the 
discrete phase-space points around the z (vertical) axis. 


This identification enables to combine sampling of the 
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initial state in the discrete phase-space and continuous 
simulation in the SU{N) phase space [27]. Further, con¬ 
tinuous phase-space methods enable a systematic expan¬ 
sion in the quantum noise parameter [TTl [12] and exten¬ 
sion to the dissipative models m- However, in order to 
use the phase-space methods one needs to find a differen¬ 
tial representation for the operators on the local Hilbert 
space. This is in general a difficult task, which was ac¬ 
complished in for any SU(N) group only in the s = 1 (P 
representation) case [33], and for any s only in the two 
dimensional case m (i.e. for the SU(2) group). Once 
the equations of motion in the phase space are obtained 
the main obstacles for their efficient simulation are di¬ 
vergence of stochastic trajectories and negativity of the 
diffusion kernel. The later can be avoided by using the 
positive-P representation (positive-P function), whereas 
the former can be to some extend circumvented by ex¬ 
ploiting a particular gauge freedom [34] . 


III. DISCRETE POSITIVE-P DISTRIBUTION, 
DIVERGING TRAJECTORIES AND NOISE 
REDUCTION 


In the rest of the paper we discuss the applicability 
of the positive-P function on discrete phase spaces and 
propose a discrete phase-space method to treat the prob¬ 
lem of diverging trajectories. The normalised positive-P 
function is usually defined through the off-diagonal co¬ 
herent state expansion of the density matrix [33] 

p = (is) 


where |U) := A(U)|0) denotes a normalized coherent 
state. We can write above equation (13) in terms of the 
quantisation kernel 0 


as 


P = 



d/i(u)d/i(o')^^"H^. u')A^-"^(u, n'), 


with the kernel 




tr (A(fi)A(*)At(ll'))' 


(14) 


(15) 


By choosing 5 = 1 in equation ([id] ) we obtain the den¬ 
sity matrix expansion given in ( |13[ ) [42] . The positive-P 
function P^+^^(U, O') always exists and can be expressed 
in terms of the discrete phase-space values as 


n') = (le) 


where p(+^)(Hc,/3) = tr (^pA^"*"^^) and S{fl) is the delta 

function for the measure /i(0). Above equation (16) 
means that the discrete P distribution can be inter¬ 
preted as a continuous positive-P distribution. Hence, 
the discrete-P representation can be used to sample the 
initial condition in simulations of the positive P function. 


Importance of the initial condition sampling and a com¬ 
parison of continuous and discrete distributions was dis¬ 
cussed in [27]. However, the derivation in [27] is valid 
only for Wigner functions of spin 1/2 systems and is 
based on a product-probability assumption, which results 
in a truncated Wigner-type approximation not permit¬ 
ting a systematic expansion or application to open quan¬ 
tum systems. The approach presented here is valid for 
any V, any s-ordered quasiprobability distribution and 
the generalised positive-P function, enables a systematic 
expansion in the noise terms, and is applicable to open 
quantum systems. 

Using the phase-space correspondence one can express 
the action of the elements of the su{N) algebra on the 
density matrix by using only first derivatives with re¬ 
spect to the phase-space variables [33]. Hence, any evo¬ 
lution equation for the density matrix which contains at 
most linear elements in the su{N) generators can be ex¬ 
pressed as an ordinary differential equation on the ex¬ 
tended phase space. On the other hand, any evolution 
equation which contains at most quadratic elements in 
the su{N) generators can be expressed as a stochastic dif¬ 
ferential equation on the extended phase space. Stochas¬ 
tic terms arise from the interaction between different sites 
of a many-body system or dissipation and make the sim¬ 
ulation of evolution equations inefficient for longer times. 
One procedure to reduce the noise is to add a gauge 
or to experiment with different decompositions of the 
diffusion matrix [34]. Another possibility is to enlarge 
the local Hilbert space with dimension N by including k 
nearest neighbours. Inside the block the evolution in the 
SU(N^) generalised phase space is described by ordinary 
differential equations, whereas the interaction between 
the blocks is still treated stochastically. The decomposi¬ 
tion into larger blocks is still exact but reduces the num¬ 
ber of stochastic terms at the expense of enlarging the 
local phase space [43] . 

Reduction of noise does not necessary solve the prob¬ 
lem of diverging stochastic trajectories, which can be 
avoided by using discrete distributions as follows. When¬ 
ever the phase-space variable becomes too large we ex¬ 
pand the corresponding kernel in terms of an equivalent 
discrete distribution. Then we randomly choose (accord¬ 
ing to the obtained discrete distribution) one of the dis¬ 
crete phase-point operators (kernels) and continue the 
simulation with the chosen kernel (for details see Ap¬ 
pendix]^ and Appendix [C]). For a large number of sim¬ 
ulated trajectories this procedure converges to the exact 
result. 


IV. EXAMPLE: DISSIPATIVE LONG-RANGE 
TRANSVERSE ISING CHAIN 


As an example of the discrete phase space projection 
method described in the previous Section we consider 
a long-range dissipative Ising chain in a transverse field 
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with the Hamiltonian 




J{a,n) 


■k\^ 






a: 




(17) 


boundary Lindblad operators Li = L 2 = 

Ls = ,/^cr+, L 4 = A/^cr+, and bulk dephasing 
j = 1, 2,... n. We are using the Kac nor¬ 
malization J{a^n) = The time evolution 

of the density matrix is given by the Lindblad equation 

H 1 

-p=-i[H,p] + J2L,pLl--{LlL„p} (18) 

At=l 

and can be expressed as a partial differential equation 
on the extended phase space [ 2 ] (see Appendix]^. This 


problem can not be solved exactly by any known method 
to treat open quantum systems [35H39] and it is hard to 
simulate with existing numerical methods [40]. Usually 
we are interested in the steady state of open quantum 
systems. In order to estimate the long-time behaviour 
we consider the interaction semiclassically, whereas the 
dissipation is treated fully quantum mechanically. For 
small system sizes (n = 5) we compare the results ob¬ 
tained by the proposed method with exact results and 
observe good agreement for short and long times. It is 
surprising that semiclassical approximation of the inter¬ 
action provides an accurate estimation of the long-time 
behaviour even in the strongly interacting regime (see 
Figure [^. We also compute the long-time averages for 
a larger system with n = 20. Parameters used in pre¬ 
sented simulations are: (a = 1.5, h = l, 71 = 0 . 2 , 72 = 
0.02, 73 = 0.1, 74 = 0.05, 7 d = 0.001. 


V. CONCLUSIONS AND DISCUSSION 

We have shown that the phase-point operators defin¬ 
ing discrete quasiprobability distributions are equivalent 
to the continuous quantisation kernel evaluated at spa- 
cial points of the phase space. This identification justi¬ 
fies use of continuous phase-space methods on discrete 
phase spaces. Further, we demonstrated that discrete 
phase-space sampling can be used to treat the problem 
of diverging trajectories. We applied the procedure on a 
dissipative transverse Ising chain and calculated the long¬ 
time expectation values of local observables up to system 
size n = 20. We found that long-time behaviour of local 
observables is well described by a semiclasical approxi¬ 
mation of the interaction even in the strongly interacting 
regime. It is an interesting question if this generalises to 
other models. For example to the integrable XXZ chain. 


where the steady state is described by a quasilocal op¬ 
erator [41]. The discrete phase space can be regarded 
as a special basis in the space of operators. Thus, the 
established connection between discrete and continuous 
phase space may enable a construction of mixed Hilbert- 
space-phase-space methods for simulation of many-body 
quantum systems. From a more theoretical point of view 
it would be interesting to explore its implications on the 
problem of mutually unbiased bases in non-prime number 
dimensions. 
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Appendix A: Derivation of the Fokker-Planck equation for the positive-P function 

In this Appendix we derive a partial-differential equation of the problem considered in paper. The Hamiltonian is 
given by 




fci 2 b' - 




(Al) 


the boundary Lindblad operators are Li = L 2 = 5 ^4 ~ bulk dephasing 

Lindblad operators are Ej +4 = j = 1,.. .n. We are using the Kac normalization J{a^n) = 

Time evolution of the density matrix is determined by the Lindblad equation 

d 1 

-p = -i[H, p] + ^ L^pLl - -{LlL^, p}. 

fl=l 


(A 2 ) 


In order to arrive at the phase space formulation of the equation (A 2 ) we follow [33] • First we define unnormalized 
SU ( 2 ) coherent states for j — th spin as 


j) — |0) + V’i |1) • 


A coherent state for the complete system is then 


= IIV’i) ® H 2 ) ® • 


0 - 


(A3) 


(A4) 


The action of pauli operators af, a = x, y, z on a coherent state can be expressed in terms of partial derivatives as 


t) = + (1 - \\t) 


(A5) 


IIU (“V’i + (1 + iU ’ 

IIU = (1 - IIU • 

We can expand the density matrix in terms of many-body coherent states as 

JM 


(A 6 ) 







where the bar • denotes complex conjugation. In order to derive the partial differential equation for the density 
matrix we have to specify the action of the Pauli operators on the qunatization kernel A('0, 0) 


~ (A7) 

ctJA(;^,^) = i + (1 + A(3A,^), 

= i + (1 + Ht’l), 

A{±,l)(j] = - A{±,f). 

Applying above equations to the commutator of the Hamiltonian with the quantisation kernel A (in the following 
we use for brevity A instead of A(3!^,^) and J instead of J{a^n)) we find 




i E 01 . ((1 - - (1 - ((1 - - (1 




■ V^/c^/c 


(A8) 


-E 


j 


fci 2 b- - fch 


((1 - V’i)(l - i’l)d4’,d^u - (1 - bj)(l - <t>l)d4>jd4>u) A 


) A. 

J -1 


Similarly we compute the action of the dissipators on the quantization kernel 


’^} = 
a|Aa| - A = 


3 + '4^j4^j 

2 H" ‘2ipj(l)j 

1 + Sipjipj 

2 “h 2 ' 0 j 0 j 
—2 “h 2'0j0j 

1 + 


(jpjdrijj. + (pjdfj).) + drijj.dfj). j A, 




(A9) 


In order to further simplify the notation we define a vector of phase space variables z 0 i, ... 0^) and 

the partial derivatives dj = 9^^., 9^+^ = 90^.. Using the equations (A8) and (A9) we obtain the following equation for 
the density matrix 


^P= [ dfi{f)dii{f)CAA, (AlO) 

JM 

2n ^ 2n 

£a = ^^Ajdj + - Dj^kdjdk: 
j=l j,k=l 
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where the non-zero elements of the drift vector A and the diffusion matrix D are given by 
3 + / I 1 + 


=^3,1 (V’i+ 72 :' 


^n+jf —( Tl 


2 “h ‘ 2 ipj(l)j 


— 2ih^pj — i 

/c=l; k^j 

3 + 


2 2 ipj(l)j 

j f i^k 0i+l 


^j,n 


3 + 1 + 3'0j0j 


2 -|- 2 'iljj(f)j 


■ 


—2 -|- 2 ipj(l)j 


1 -h l^l)j<pj 


(All) 


, V i + ^/’fcbfc 

1 + 3'ipj(j)j 


(1 + ) ; 


■ \2h(j)j + i 


'i + 72 
J 


2 “h 2ipj(l)j ^ 2 2'0j0j 


( 73 


3 + 


7 + 74 


1 + ^'ipjcjyj 


2 “h 2'ipj(l)j 2 “h 2ipj(l)j 


/ + 7d 


—2 “h 2 ipj(l)j 

1 -h 


i^k + 0/c 

2|j - Vl + 


(1+0^) , 


/c=l; k^j 

Dj+n,j =Dj,j+n = (5i,l(7l + l2'lpj4>])Sj,ni'y3 + +i'yD'lpj4>j, 

■^i.i+1 ~^3+i,j ~ ~or7~TTa ~ ~ 7 + 1 )’ 


2|j - fch 


a 


'n+j,n+j+l —Dn-\-j-\-l,n+j 


Ji 


2\j-k\^ 




By partial integration of the equation (AlO) and assuming that boundary terms vanish we obtain the following partial 
differential equation for the positive-P function 


d 

dt 


P{z) = 


^ 2n . 2n \ 

-E77 + 2 E djdkDj,A P{z). 

^ i=i i,fe=i / 


(A 12 ) 


This equation is a Fokker-Planck equation which can be simulated by rewriting it to a stochastic differential equation 
(see e.g. [33]). 


Appendix B: Discrete phase-space solution for divergent trajectories 


In this section a simple example of discrete phase-space projection is discussed. We will demonstrate that the 
diverging trajectories are a consequence of a particular choice of the phase space and could in principle be avoided by 
using a redundant representation of the phase space. We consider the SU(2) case (i.e. a spin 1/2). In this case the 
phase space pertaining to the positive-P function is C^, which can be conveniently represented by using stereographic 
projection to two spheres. To each point in the phase space corresponds a kernel A{z,v) = \\z) (f;|| / {^z). Since 
this is an operator in the Hilbert space it can be expanded using the discrete phase space points as A{z, v) = 
Sa /3 7 ( 5 Pa,/ 3 , 7 ,^A( 2 :Q,^/ 3 , 2 :^,^), with z^^f^ defined in Section [ll| and Pa,p,-f,6 positive real numbers summing to one. If at 
some point in the stochastic evolution \zj{t)\ or \vj{t)\ crosses a threshold z^ax we calculate the discrete probability 
distribution ^ ^ corresponding to the phase-space variables Zj{t)^ time t and randomly choose one discrete 

phase space point (za^p^Vj^s) according to the distribution Pa,p,j,5- We interpret the chosen kernel (phase space point 
{zj{t),Vj{t)) as one corresponding to a continuous-phase space distribution and continue the evolution 

from this point. Phase-space variables corresponding to particles at sites /c 7^ j do not change. 

We found numerically that in all cases Pa,^,j,6 are positive and sum to one. In Figure we show an example 
trajectory corresponding to the variable zi{t) of the first spin. Since we fixed the highest weight representation the 
condition \z{t)\ < Zmax means that the spin in the z direction cannot be smaller than cos{0max/‘^)' If at some time 
the trajectory of the phase-space variable crosses the critical circle on the sphere we project the phase-space point at 
the crossing to one point on the sphere corresponding to its discrete quasi-probability distribution. From this picture 
it is clear that the divergence of the trajectory is related to the inability of the chosen kernel to represent all states in 
the Hilbert space by coherent states corresponding to the phase-space region with |2)| < Zmax- 

In Figure|^we show the average number of projections (jumps to the discrete distribution) per trajectory for system 
sizes n = 5, 20. 
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Figure 3: Demonstration of the algorithm on a particular trajectory. The red and blue points denote the stochastic trajectory 
before and after the jump, respectively. The brightness of the color denotes how far in the past is the trajectory (the initial 
point and the point after the jump are white). The dark red points are close to the Zmax = 10a/ 2, when the trajectory crosses 
this value it is projected to the blue discrete phase-space point. After that the stochastic evolution continues. 



Figure 4: Average number of discrete projections per trajectory until time t. The black, blue dots correspond to n = 5, 20. 
The full line serves as a guide to the eye and is proportional to . 


Appendix C: Calculation of local expectation values 

Once we have an ensemble of stochastic trajectories we calculate the expectation values of local observables (e.g. 
on site j) by stochastic averages 


(Oj) = tr{pOj) = {Oi'ipjAj)) 

where (•)stoch denotes an average over stochastic trajectories and 

In the spin 1/2 case the local observables are given by 


I'/’jV’i - 




1 — 7/0 


1 + 07 / ’ ’ 1 + 00 ’ ’ 1 + 00 

= — 1 we perform discrete projections also when 


(Cl) 

(C2) 


(C3) 


Since these functions have poles at 00 
simulations we take e = 0.1). 


1| < e (in the 









